*************************************************************
*** Analyze the Monte Carlo experiments (SEM data-generating process) for the WWW JoP project
***
*** Created: 9-19-14
*** Modified: 9-5-19
***
*************************************************************

*** Open the data containing the SEM simulations
use "Experiments\SEM\SEM.dta", clear

preserve
	duplicates drop lambda, force
	gen id = _n
	local end = _N
	sort lambda
	tempfile id
	save `id', replace 
restore

sort lambda
merge lambda using `id'
drop _merge

tempname sem
postfile `sem' str20 model lambda rec_beta rec_lambda avg_beta avg_lambda /*
*/	rec_teffect rec_teffect_o1 rec_teffect_o2 rec_teffect_o3 rec_deffect rec_deffect_o0 rec_deffect_o1 rec_deffect_o2 rec_deffect_o3 rec_ieffect rec_ieffect_o1 rec_ieffect_o2 rec_ieffect_o3 rec_feffect rec_heffect /*
*/	avg_teffect avg_teffect_o1 avg_teffect_o2 avg_teffect_o3 avg_deffect avg_deffect_o0 avg_deffect_o1 avg_deffect_o2 avg_deffect_o3 avg_ieffect avg_ieffect_o1 avg_ieffect_o2 avg_ieffect_o3 avg_feffect avg_heffect /*
*/	true_teffect true_deffect true_ieffect true_feffect true_heffect true_deffect_o0 true_teffect_o1 true_deffect_o1 true_ieffect_o1 true_teffect_o2 true_deffect_o2 true_ieffect_o2 true_teffect_o3 true_deffect_o3 true_ieffect_o3 /*
*/	using "Experiments\Analyze\SEM\Analyze SEM Monte Carlo Experiments.dta", replace 


*** SEM model
preserve 
	gen sig_b_x = cond((b_x_sem - 1.96*se_x_sem < b1) & (b_x_sem + 1.96*se_x_sem > b1), 1, 0)
	gen sig_lambda = cond((lambda_sem - 1.96*se_lambda_sem < lambda) & (lambda_sem + 1.96*se_lambda_sem > lambda), 1, 0)
	gen sig_teffect = cond((b_x_sem - 1.96*se_x_sem < true_teffect) & (b_x_sem + 1.96*se_x_sem > true_teffect), 1, 0)
	gen sig_deffect = cond((b_x_sem - 1.96*se_x_sem < true_deffect) & (b_x_sem + 1.96*se_x_sem > true_deffect), 1, 0)
	gen sig_deffect_o0 = cond((deffect_o0_sem - 1.96*deffect_o0_se_sem < true_deffect) & (deffect_o0_sem + 1.96*deffect_o0_se_sem > true_deffect), 1, 0)
	
	qui foreach n of numlist 1(1)`end' {
		qui sum lambda if id == `n'
		local r = round(r(mean), .001)

		foreach v of varlist true_teffect - true_ieffect_o3 {
			qui sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}

		foreach v in teffect deffect deffect_o0 b_x lambda {
			sum `v'_sem if id == `n', meanonly
			local avg_`v' = r(mean)
			
			sum sig_`v' if id == `n', meanonly
			local rec_`v' = 100 * r(mean)
		}
		
		nois di _newline(2) "******************************************************"
		nois di "Lambda = `r'; Beta_x = 1"
		nois di "% of cases where SEM Model recovers the true beta_x (reject H_o: beta_x = 1) = " `rec_b_x'
		nois di "% of cases where SEM Model recovers the true lambda (reject H_o: lambda = `r') = " `rec_lambda'
		nois di _newline(2)
		
		post `sem' ("SEM") (`r') (`rec_b_x') (`rec_lambda') (`avg_b_x') (`avg_lambda') /*
		*/	(`rec_teffect') (.) (.) (.) (`rec_deffect') (`rec_deffect_o0') (.) (.) (.) (.) (.) (.) (.) (.) (.) /*
		*/	(`avg_teffect') (.) (.) (.) (`avg_deffect') (`avg_deffect_o0') (.) (.) (.) (.) (.) (.) (.) (.) (.) /*
		*/	(`true_teffect') (`true_deffect') (`true_ieffect') (`true_feffect') (`true_heffect') (`true_deffect_o0') (`true_teffect_o1') (`true_deffect_o1') (`true_ieffect_o1') (`true_teffect_o2') (`true_deffect_o2') (`true_ieffect_o2') (`true_teffect_o3') (`true_deffect_o3') (`true_ieffect_o3')		
	}
restore


*** SAR model
preserve
	gen sig_b = cond((b_x_sar - 1.96*se_x_sar < b1) & (b_x_sar + 1.96*se_x_sar > b1), 1, 0)

	qui foreach n of numlist 1(1)`end' {
		sum lambda if id == `n'
		local r = round(r(mean), .001)

		foreach v of varlist true_teffect - true_ieffect_o3 {
			sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}
		
		sum b_x_sar if id == `n', meanonly
		local avg_beta = r(mean)

		sum deffect_o0_sar if id == `n', meanonly
		local avg_deffect_o0 = r(mean)
		
		tempvar sig_deffect_o0
		gen `sig_deffect_o0' = cond((deffect_o0_sar - 1.96*deffect_o0_se_sar < true_deffect_o0) & (deffect_o0_sar + 1.96*deffect_o0_se_sar > true_deffect_o0), 1, 0)
		sum `sig_deffect_o0' if id == `n'
		local rec_deffect_o0 = 100 * r(mean)
		
		foreach e in f h {
			sum `e'effect_sar if id == `n', meanonly
			local avg_`e'effect = r(mean)
			
			tempvar sig_`e'effect
			gen `sig_`e'effect' = cond((`e'effect_sar - 1.96*`e'effect_se_sar < true_`e'effect) & (`e'effect_sar + 1.96*`e'effect_se_sar > true_`e'effect), 1, 0)
			sum `sig_`e'effect' if id == `n', meanonly
			local rec_`e'effect = 100 * r(mean)
		}
		
		foreach v in effect effect_o1 effect_o2 effect_o3 {
			foreach e in t d i {
				sum `e'`v'_sar if id == `n', meanonly
				local avg_`e'`v' = r(mean)
			
				tempvar sig_`e'`v' 
				gen `sig_`e'`v'' = cond((`e'`v'_sar - 1.96*`e'`v'_se_sar < true_`e'`v') & (`e'`v'_sar + 1.96*`e'`v'_se_sar > true_`e'`v'), 1, 0)
				sum `sig_`e'`v'' if id == `n', meanonly
				local rec_`e'`v' = 100 * r(mean)
			}		
		}
		
		foreach v in b {
			sum sig_`v' if id == `n'
			local rec_`v' = 100 * r(mean)
		}
		
		nois di _newline(2) "******************************************************"
		nois di "Lambda = `r'; Beta_x = 1"
		nois di "% of cases where SAR Model recovers the average total effect of x (beta_x = 1) = " `rec_teffect'
		nois di "% of cases where SAR Model recovers the average first-order total effect of x (beta_x = 1) = " `rec_teffect_o1'
		nois di "% of cases where SAR Model recovers the average second-order total effect of x (beta_x = 1) = " `rec_teffect_o2'
		nois di "% of cases where SAR Model recovers the average third-order total effect of x (beta_x = 1) = " `rec_teffect_o3'		
		nois di _newline(1)

		nois di "% of cases where SAR Model recovers the average direct effect of x (beta_x = 1) = " `rec_deffect'
		nois di "% of cases where SAR Model recovers the average zero-order direct effect of x (beta_x = 1) = " `rec_deffect_o0'
		nois di "% of cases where SAR Model recovers the average second-order direct effect of x (beta_x = 1) = " `rec_deffect_o2'
		nois di "% of cases where SAR Model recovers the average third-order direct effect of x (beta_x = 1) = " `rec_deffect_o3'		
		nois di _newline(1)
		
		nois di "% of cases where SAR Model recovers the average indirect effect of x (beta_x = 1) = " `rec_ieffect'
		nois di "% of cases where SAR Model recovers the average first-order indirect effect of x (beta_x = 1) = " `rec_ieffect_o1'
		nois di "% of cases where SAR Model recovers the average second-order indirect effect of x (beta_x = 1) = " `rec_ieffect_o2'
		nois di "% of cases where SAR Model recovers the average third-order indirect effect of x (beta_x = 1) = " `rec_ieffect_o3'		
		nois di _newline(1)
		
		nois di "% of cases where SAR Model recovers the true higher-order effect of x (beta_x = 1) = " `rec_heffect'
		nois di "% of cases where SAR Model recovers the true feedback effect of x (beta_x = 1) = " `rec_feffect'
		nois di "% of cases where SAR Model falsely rejects the true null of no higher-order effects = " 100 - `rec_heffect'
		nois di "% of cases where SAR Model falsely rejects the true null of no feedback effects = " 100 - `rec_feffect'
		nois di _newline(2)
		
		post `sem' ("SAR") (`r') (`rec_b') (.) (`avg_beta') (.) /*
		*/	(`rec_teffect') (`rec_teffect_o1') (`rec_teffect_o2') (`rec_teffect_o3') (`rec_deffect') (`rec_deffect_o0') (`rec_deffect_o1') (`rec_deffect_o2') (`rec_deffect_o3') (`rec_ieffect') (`rec_ieffect_o1') (`rec_ieffect_o2') (`rec_ieffect_o3') (`rec_feffect') (`rec_heffect') /*
		*/	(`avg_teffect') (`avg_teffect_o1') (`avg_teffect_o2') (`avg_teffect_o3') (`avg_deffect') (`avg_deffect_o0') (`avg_deffect_o1') (`avg_deffect_o2') (`avg_deffect_o3') (`avg_ieffect') (`avg_ieffect_o1') (`avg_ieffect_o2') (`avg_ieffect_o3') (`avg_feffect') (`avg_heffect') /*
		*/	(`true_teffect') (`true_deffect') (`true_ieffect') (`true_feffect') (`true_heffect') (`true_deffect_o0') (`true_teffect_o1') (`true_deffect_o1') (`true_ieffect_o1') (`true_teffect_o2') (`true_deffect_o2') (`true_ieffect_o2') (`true_teffect_o3') (`true_deffect_o3') (`true_ieffect_o3')		
	}
restore


*** SLX model 1: y = BX + Theta W X
preserve
	local m = 1

	tempvar sig_b_x sig_deffect_o0
	gen `sig_b_x' = cond((b_x_m`m' - 1.96*se_x_m`m' < b1) & (b_x_m`m' + 1.96*se_x_m`m' > b1), 1, 0)
	gen `sig_deffect_o0' = cond((deffect_o0_m`m' - 1.96*deffect_o0_se_m`m' < true_deffect_o0) & (deffect_o0_m`m' + 1.96*deffect_o0_se_m`m' > true_deffect_o0), 1, 0)
	
	foreach e in t d i {
		tempvar sig_`e'effect 
		gen `sig_`e'effect' = cond((`e'effect_m`m' - 1.96*`e'effect_se_m`m' < true_`e'effect) & (`e'effect_m`m' + 1.96*`e'effect_se_m`m' > true_`e'effect), 1, 0)

		foreach o in 1 {
			tempvar sig_`e'effect_o`o'
			gen `sig_`e'effect_o`o'' = cond((`e'effect_o`o'_m`m' - 1.96*`e'effect_o`o'_se_m`m' < true_`e'effect_o`o') & (`e'effect_o`o'_m`m' + 1.96*`e'effect_o`o'_se_m`m' > true_`e'effect_o`o'), 1, 0)
		}
	}
	
	qui foreach n of numlist 1(1)`end' {
		sum lambda if id == `n', meanonly
		local r = round(r(mean), .001)

		foreach v of varlist true_teffect - true_ieffect_o3 {
			sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}
		
		foreach v in b_x teffect teffect_o1 deffect deffect_o0 deffect_o1 ieffect ieffect_o1 {
			sum `v'_m`m' if id == `n', meanonly
			local avg_`v' = r(mean)
			
			sum `sig_`v'' if id == `n', meanonly
			local rec_`v' = 100 * r(mean)
		}
		
		nois di _newline(2) "******************************************************"
		nois di "Lambda = `r'; Beta_x = 1"
		nois di "% of cases where SLX Model (1 W) recovers the average total effect = " `rec_teffect'
		nois di "% of cases where SLX Model (1 W) recovers the average first-order total effect = " `rec_teffect_o1'
		nois di _newline(1) "% of cases where SLX Model (1 W) recovers the total average direct effect = " `rec_deffect'
		nois di "% of cases where SLX Model (1 W) recovers the average zero-order direct effect = " `rec_deffect_o0'
		nois di _newline(1) "% of cases where SLX Model (1 W) recovers the total average indirect effect = " `rec_ieffect'
		nois di "% of cases where SLX Model (1 W) recovers the average first-order indirect effect = " `rec_ieffect_o1'
		nois di _newline(2)
		
		post `sem' ("SLX1") (`r') (`rec_b_x') (.) (`avg_b_x') (.) /*
		*/	(`rec_teffect') (`rec_teffect_o1') (.) (.) (`rec_deffect') (`rec_deffect_o0') (`rec_deffect_o1') (.) (.) (`rec_ieffect') (`rec_ieffect_o1') (.) (.) (.) (.) /*
		*/	(`avg_teffect') (`avg_teffect_o1') (.) (.) (`avg_deffect') (`avg_deffect_o0') (`avg_deffect_o1') (.) (.) (`avg_ieffect') (`avg_ieffect_o1') (.) (.) (.) (.) /*
		*/	(`true_teffect') (`true_deffect') (`true_ieffect') (`true_feffect') (`true_heffect') (`true_deffect_o0') (`true_teffect_o1') (`true_deffect_o1') (`true_ieffect_o1') (`true_teffect_o2') (`true_deffect_o2') (`true_ieffect_o2') (`true_teffect_o3') (`true_deffect_o3') (`true_ieffect_o3')		
	}
restore


*** SLX model 2: y = BX + Theta_1 W X + Theta_2 W^2 X
preserve
	local m = 2

	tempvar sig_b_x sig_deffect_o0
	gen `sig_b_x' = cond((b_x_m`m' - 1.96*se_x_m`m' < b1) & (b_x_m`m' + 1.96*se_x_m`m' > b1), 1, 0)
	gen `sig_deffect_o0' = cond((deffect_o0_m`m' - 1.96*deffect_o0_se_m`m' < true_deffect_o0) & (deffect_o0_m`m' + 1.96*deffect_o0_se_m`m' > true_deffect_o0), 1, 0)
	
	foreach v in teffect deffect ieffect feffect heffect {
		tempvar sig_`v' 
		gen `sig_`v'' = cond((`v'_m`m' - 1.96*`v'_se_m`m' < true_`v') & (`v'_m`m' + 1.96*`v'_se_m`m' > true_`v'), 1, 0)
	}
	
	foreach e in t d i {
		foreach o of numlist 1(1)2 {
			tempvar sig_`e'effect_o`o'
			gen `sig_`e'effect_o`o'' = cond((`e'effect_o`o'_m`m' - 1.96*`e'effect_o`o'_se_m`m' < true_`e'effect_o`o') & (`e'effect_o`o'_m`m' + 1.96*`e'effect_o`o'_se_m`m' > true_`e'effect_o`o'), 1, 0)
		}
	}

	qui foreach n of numlist 1(1)`end' {
		qui sum lambda if id == `n'
		local r = round(r(mean), .001)

		foreach v of varlist true_teffect - true_ieffect_o3 {
			qui sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}
		
		foreach v in b_x teffect teffect_o1 teffect_o2 deffect deffect_o0 deffect_o1 deffect_o2 ieffect ieffect_o1 ieffect_o2 feffect heffect {
			sum `v'_m`m' if id == `n', meanonly
			local avg_`v' = r(mean)
			
			sum `sig_`v'' if id == `n', meanonly
			local rec_`v' = 100 * r(mean)
		}
		
		nois di _newline(2) "******************************************************"
		nois di "Labmda = `r'; Beta_x = 1"
		nois di "% of cases where SLX Model (2 W) recovers the average total effect = " `rec_teffect'
		nois di "% of cases where SLX Model (2 W) recovers the average first-order total effect = " `rec_teffect_o1'
		nois di "% of cases where SLX Model (2 W) recovers the average second-order total effect = " `rec_teffect_o2'	
		nois di _newline(1) "% of cases where SLX Model (2 W) recovers the total average direct effect = " `rec_deffect'
		nois di "% of cases where SLX Model (2 W) recovers the average zero-order direct effect = " `rec_deffect_o0'
		nois di "% of cases where SLX Model (2 W) recovers the average second-order direct effect = " `rec_deffect_o2'
		nois di _newline(1) "% of cases where SLX Model (2 W) recovers the total average indirect effect = " `rec_ieffect'
		nois di "% of cases where SLX Model (2 W) recovers the average first-order indirect effect = " `rec_ieffect_o1'
		nois di "% of cases where SLX Model (2 W) recovers the average second-order indirect effect = " `rec_ieffect_o2'
		nois di "% of cases where SLX Model (2 W) recovers the average feedback effects = " `rec_feffect'
		nois di "% of cases where SLX Model (2 W) recovers the average higher-order effects = " `rec_heffect'		
		nois di _newline(2)
		
		post `sem' ("SLX2") (`r') (`rec_b_x') (.) (`avg_b_x') (.) /*
		*/	(`rec_teffect') (`rec_teffect_o1') (`rec_teffect_o2') (.) (`rec_deffect') (`rec_deffect_o0') (`rec_deffect_o1') (`rec_deffect_o2') (.) (`rec_ieffect') (`rec_ieffect_o1') (`rec_ieffect_o2') (.) (`rec_feffect') (`rec_heffect') /*
		*/	(`avg_teffect') (`avg_teffect_o1') (`avg_teffect_o2') (.) (`avg_deffect') (`avg_deffect_o0') (`avg_deffect_o1') (`avg_deffect_o2') (.) (`avg_ieffect') (`avg_ieffect_o1') (`avg_ieffect_o2') (.) (`avg_feffect') (`avg_heffect') /*
		*/	(`true_teffect') (`true_deffect') (`true_ieffect') (`true_feffect') (`true_heffect') (`true_deffect_o0') (`true_teffect_o1') (`true_deffect_o1') (`true_ieffect_o1') (`true_teffect_o2') (`true_deffect_o2') (`true_ieffect_o2') (`true_teffect_o3') (`true_deffect_o3') (`true_ieffect_o3')		
	}
restore


*** SLX model 3: y = BX + Theta_1 W X + Theta_2 W^2 X + Theta_3 W^3 X
preserve
	local m = 4

	tempvar sig_b_x sig_deffect_o0
	gen `sig_b_x' = cond((b_x_m`m' - 1.96*se_x_m`m' < b1) & (b_x_m`m' + 1.96*se_x_m`m' > b1), 1, 0)
	gen `sig_deffect_o0' = cond((deffect_o0_m`m' - 1.96*deffect_o0_se_m`m' < true_deffect_o0) & (deffect_o0_m`m' + 1.96*deffect_o0_se_m`m' > true_deffect_o0), 1, 0)
	
	foreach v in teffect deffect ieffect feffect heffect {
		tempvar sig_`v' 
		gen `sig_`v'' = cond((`v'_m`m' - 1.96*`v'_se_m`m' < true_`v') & (`v'_m`m' + 1.96*`v'_se_m`m' > true_`v'), 1, 0)
	}
	
	foreach e in t d i {
		foreach o of numlist 1(1)3 {
			tempvar sig_`e'effect_o`o'
			gen `sig_`e'effect_o`o'' = cond((`e'effect_o`o'_m`m' - 1.96*`e'effect_o`o'_se_m`m' < true_`e'effect_o`o') & (`e'effect_o`o'_m`m' + 1.96*`e'effect_o`o'_se_m`m' > true_`e'effect_o`o'), 1, 0)
		}
	}
	
	qui foreach n of numlist 1(1)`end' {
		qui sum lambda if id == `n'
		local r = round(r(mean), .001)

		foreach v of varlist true_teffect - true_ieffect_o3 {
			qui sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}
		
		foreach v in b_x teffect teffect_o1 teffect_o2 teffect_o3 deffect deffect_o0 deffect_o1 deffect_o2 deffect_o3 ieffect ieffect_o1 ieffect_o2 ieffect_o3 feffect heffect {
			sum `v'_m4 if id == `n', meanonly
			local avg_`v' = r(mean)
			
			sum `sig_`v'' if id == `n', meanonly
			local rec_`v' = 100 * r(mean)
		
		}	
		
		nois di _newline(2) "******************************************************"
		nois di "Lambda = `r'; Beta_x = 1"
		nois di "% of cases where SLX Model (3 W) recovers the average total effect = " `rec_teffect'
		nois di "% of cases where SLX Model (3 W) recovers the average first-order total effect = " `rec_teffect_o1'
		nois di "% of cases where SLX Model (3 W) recovers the average second-order total effect = " `rec_teffect_o2'
		nois di "% of cases where SLX Model (3 W) recovers the average third-order total effect = " `rec_teffect_o3'
		nois di _newline(1) "% of cases where SLX Model (3 W) recovers the total average direct effect = " `rec_deffect'
		nois di "% of cases where SLX Model (3 W) recovers the average zero-order direct effect = " `rec_deffect_o0'
		nois di "% of cases where SLX Model (3 W) recovers the average second-order direct effect = " `rec_deffect_o2'
		nois di "% of cases where SLX Model (3 W) recovers the average third-order direct effect = " `rec_deffect_o3'
		nois di _newline(1) "% of cases where SLX Model (3 W) recovers the total average indirect effect = " `rec_ieffect'
		nois di "% of cases where SLX Model (3 W) recovers the average first-order indirect effect = " `rec_ieffect_o1'
		nois di "% of cases where SLX Model (3 W) recovers the average second-order indirect effect = " `rec_ieffect_o2'
		nois di "% of cases where SLX Model (3 W) recovers the average third-order indirect effect = " `rec_ieffect_o3'
		nois di "% of cases where SLX Model (3 W) recovers the average feedback effects = " `rec_feffect'
		nois di "% of cases where SLX Model (3 W) recovers the average higher-order effects = " `rec_heffect'		
		nois di _newline(2)
		
		post `sem' ("SLX3") (`r') (`rec_b_x') (.) (`avg_b_x') (.) /*
		*/	(`rec_teffect') (`rec_teffect_o1') (`rec_teffect_o2') (`rec_teffect_o3') (`rec_deffect') (`rec_deffect_o0') (`rec_deffect_o1') (`rec_deffect_o2') (`rec_deffect_o3') (`rec_ieffect') (`rec_ieffect_o1') (`rec_ieffect_o2') (`rec_ieffect_o3') (`rec_feffect') (`rec_heffect') /*
		*/	(`avg_teffect') (`avg_teffect_o1') (`avg_teffect_o2') (`avg_teffect_o3') (`avg_deffect') (`avg_deffect_o0') (`avg_deffect_o1') (`avg_deffect_o2') (`avg_deffect_o3') (`avg_ieffect') (`avg_ieffect_o1') (`avg_ieffect_o2') (`avg_ieffect_o3') (`avg_feffect') (`avg_heffect') /*
		*/	(`true_teffect') (`true_deffect') (`true_ieffect') (`true_feffect') (`true_heffect') (`true_deffect_o0') (`true_teffect_o1') (`true_deffect_o1') (`true_ieffect_o1') (`true_teffect_o2') (`true_deffect_o2') (`true_ieffect_o2') (`true_teffect_o3') (`true_deffect_o3') (`true_ieffect_o3')		
	}
restore


postclose `sem'

*************************************************************************
************************* Analyze ***************************************
*************************************************************************
use "Experiments\Analyze\SEM\Analyze SEM Monte Carlo Experiments.dta", clear

***************************** True Effects (SEM) *************************
list true_teffect - true_heffect if model == "SEM"

***************************** Recovery Rates *****************************
*** Recovery rates of SAR: y = XB + pWy
list lambda rec_teffect - rec_heffect if model == "SAR"

*** Recovery rates of SLX 1: y = XB + Theta W X
list lambda rec_teffect - rec_heffect if model == "SLX1"

*** Recovery rates of SLX 2: y = XB + Theta_1 W X + Theta_2 W^2 X
list lambda rec_teffect - rec_heffect if model == "SLX2"

*** Recovery rates of SLX 3: y = XB + Theta_1 W X + Theta_2 W^2 X + Theta_3 W^3 X
list lambda rec_teffect - rec_heffect if model == "SLX3"

***************************** Average Effects *****************************

*** Average effects of SAR: y = XB + pWy
list lambda avg_teffect - avg_heffect if model == "SAR"

*** Average effects of SLX 1: y = XB + Theta W X
list lambda avg_teffect - avg_heffect if model == "SLX1"

*** Average effects of SLX 2: y = XB + Theta_1 W X + Theta_2 W^2 X
list lambda avg_teffect - avg_heffect if model == "SLX2"

*** Average effects of SLX 3: y = XB + Theta_1 W X + Theta_2 W^2 X + Theta_3 W^3 X
list lambda avg_teffect - avg_heffect if model == "SLX3"


***************************** Higher-Order Effects *****************************

gen reject_feffect = 100 - rec_feffect
gen reject_heffect = 100 - rec_heffect

*** False rejection rates for feedback effects:
list model lambda reject_feffect avg_feffect if inlist(model, "SAR", "SLX2", "SLX3")

*** False rejection rates for higher-order effects:
list model lambda reject_heffect avg_heffect if inlist(model, "SAR", "SLX2", "SLX3")
